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Abstract 

We describe the physical model, numerical algorithms, and software structure of WRF-Fire. 
WRF-Fire consists of a fire-spread model, implemented by the level-set method, coupled 
with the Weather Research and Forecasting model. In every time step, the fire model inputs 
the surface wind, which drives the fire, and outputs the heat flux from the fire into the at- 
mosphere, which in turn influences the atmosphere. The level-set method allows submesh 
representation of the burning region and flexible implementation of various kinds of igni- 
tion. WRF-Fire is distributed as a part of WRF and it uses the WRF parallel infrastructure 
for parallel computing. i 



1 Introduction 



Wildland fires impact the lives of millions of people and cause major damage every year 
worldwide, yet they are a natural part of the cycle of nature. Better tools for modeling wild- 
land fire behavior are important for managing fire suppression, planning controlled burns to 
reduce the fuels, as well as to help assess fire danger. Fire models range from tools based 
on Rothermel ( |1972| ) fire spread rate formulas, such as BehavePlus ( Andrews") 20071 and 
FARSITE ( Finney I 1998[ l, suitable for operational forecasting, to sophisticated 3-D com- 
putational fluid dynamics and combustion simulations suitable for research and reanalysis, 
such as FIRETEC ( Linn etaL] [20021 ) and WFDS ( Mell et al.|[2007] l. BehavePlus, the PC- 
based successor of the calculator-based BEHAVE, determines the fire spread rate at a single 
point from fuel and environmental data; FARSITE uses the fire spread rate to provide a 2-D 
simulation on a PC; while FIRETEC and WFDS require a parallel supercomputer and run 
much slower than real time. 

Wildland fire is a complicated multiscale process, from the flame reaction zone on milime- 
ter scale to the synoptic weather scale of hundreds of kilometers. Since direct numerical 
simulation of wildland fire is computationally intractable and detailed data are not avail- 
able anyway, compromises in the choice of processes to be modeled, approximations, and 
parametrizations are essential. Fortunately, a practically important range of wildland fire 
behavior can be captured by the coupling of a mesoscale weather model with a simple 2-D 
fire spread model ( [Clark et al.[ 1996a|b ). Weather has a major influence on wildfire be- 
havior; in particular, wind plays a dominant role in the fire spread. Conversely, the fire 
influences the atmosphere through the heat and vapor fluxes from burning hydrocarbons 
and evaporation of fuel moisture. Fire heat output has a major effect on the atmosphere; the 
buoyancy created by the heat from the fire can cause tornadic strength winds, and the air 
motion and moisture from the fire can affect the atmosphere also away from the fire. It is 
well known that a large fire "creates its own weather". The correct wildland fire shape and 



progress result from the two-way interaction between the fire and the atmosphere (Clark 
etall |1996a'b! '2004'; 'Coin'. "2005^. 



WRF-Fire (Mandel et al. , ,2009^ combines the Weather Research and Forecasting Model 
(WRF) with the ARW dynamical core ( Skamarock et al. 2008 1 with a semi-empirical fire 
spread model. It is intended to be faster than real time in order to deliver a prediction. 

WRF-Fire has grown out of the NCAR's CAWFE code ( [Clark et al.| [I996a|b[ [2004} 



Coen 2005 1. CAWFE consists of the Clark-Hall mesoscale atmospheric model, coupled 



with a tracer-based fire spread model. Although the Clark-Hall model has many good prop- 
erties, it is a legacy serial code, not supported, and difficult to modify or use with real data, 
while WRF is a parallel supported community code routinely used for real runs. See CoenJ 
and Patton ( 2010, ) for a further discussion of their relative merits in the wildland fire appli- 
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cation. WRF-Fire was started by Patton and Coen ( 2004 1, who proposed a combination of 
WRF with the tracer-based model from CAWFE, formulated a road map, and made the im- 
portant observation that the innermost domain of the weather code, which interacts directly 
with the fire model, needs to run in the Large Eddy Simulation (LES) mode. Patton ported 
the Fortran 77-based fire module to Fortran 90 and developed the initial serial coupled 
code. However, instead of using the existing tracer-based CAWFE code, the fire module in 
WRF-Fire was developed based on the level-set method (Osher and Fedkiw 2003 1. One 
of the reasons was that the representation of the fire region by the level-set function was 
thought to be more flexible than the representation of the burning region in CAWFE by four 
tracers in each cell of the fire mesh. In particular, the level-set function can be manipulated 
more easily than tracers for the purpose of data assimilation. Insertion of the heat fluxes, 
while fundamentally the same as in CAWFE, had to be redone for WRF variables already 
in Patton's initial code. Thus, only the code for the calculation of the fire spread rate and 
the heat fluxes remained from CAWFE. While WRF-Fire takes advantage of the experience 
accumulated with CAWFE, WRF is quite different from the Clark-Hall atmospheric model 
and the fireline propagation algorithm is also different. Thus, it needs to be demonstrated 
that WRF-Fire can deliver similar results as CAWFE, and WRF-Fire needs to be validated 
against real fires (Sect. 10 1. 

The level-set method was used for a surface fire spread model in Mallet et al. ( 2009[ ). Fil- 



ippi et al. ( 2009| ) coupled the atmospheric model Meso-nh with fire propagation by tracers. 



Tiger ( Mazzoleni and Giannino[ 2010) uses a 2-D combusion model based on reaction- 
convection-diffusion equations and a convection model to emulate the effect of the fire on 
the wind. FIRESTAR (jMorvan and Dupuy 2004) is a physically accurate wildland fire 
model in two dimensions, one horizontal and one vertical. UU LES -Fire ( |Sun et aL| |2009| ) 
couples the University of Utah's Large Eddy Simulation code with the tracer-based code 
from CAWFE. See the survey bylSuUivan ( |2009 ) for a number of other models. 

WRF-Fire was briefly treated as one of the topics in Mandel et al. (2009). The pur- 
pose of this paper is to describe the fire module and the coupling with WRF in the current 
WRF-Fire code in sufficient detail, yet understandable to a reader not familiar with WRF. 
In addition, the advances since the paper Mandel et al. (2009|) was written in 2007 include 
new, practically important ignition schemes (Sect. [44] ), vertical interpolation of the wind in 
the boundary layer dependent on land-use (Sect. [6]), parallel computing (Sect. |7]), input of 
real data (Sect.|9]), and validations in progress on real fires (Sect. 10). This paper also con- 
tains reproducible descriptions of the physical model (Sect. [2]), the required WRF settings 
(Sect. [5]), and the coupling of WRF with the fire module (Sect. [6]). 

WRF-Fire is public domain software and it has been distributed as a part of the WRF 
source code at|vyrf-modeI.org| since version 3.2, released in April 2010 (|Dudhia[|2010|). The 
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released version is updated periodically and supported by NCAR. The current development 
version of WRF-Fire with the latest features and bug fixes and additional visualization tools, 
guides, and diagnostic utilities, are available directly from the developers at openwfm.orgl 
This article describes WRF-Fire as it is scheduled to be included in WRF 3.3, to be released 
in March or April 2011. WRF-Fire user's guide is available as a part of the WRF user's 



guide ( Wang et aLj 2010 1, to be updated with the release. 

New features since WRF version 3.2 include new ignition models, vertical interpolation 
of the wind from logarithmic profile, fetching high-resolution geogrid data, terrain gradient 
interpolation, and optional input of fuel map, land use map, and high-resolution topography 
in ideal runs. 



2 Physical fire model and fuels 

The physical model consists of functions specifying the fire spread rate and the heat fluxes, 
and it is essentially the same as a subset of CAWFE ( Clark et aL| [2004 Coen 2005 1. The 
spread rate calculation is in turn based on BEHAVE ( Rothermel} 1972 Andrews 2007 1. It 
is described here in more detail for the sake of reproducibility and to point out the (minor) 
differences. 



2.1 Fuel properties 

Fuel is characterized by the quantities listed in Table [T] which are given at every point of 
the fire mesh. To simplify the specification of fuel properties, fuels are given as one of 13 
Anderson ( 1982!) categories, which are preset vectors of values of the fuel properties. These 
values are specified in an input text file (namelist . fire), and they can be modified by 
the user. The user can also specify completely new, custom fuel categories. 

2.2 Fire spread rate 

The fire model is posed in the horizontal (x, y) plane the Earth surface is projected on. The 
semi-empirical approach to fire propagation used here assumes that the fire spread rate is 



given by the modified Rothermel] ( | 1 972) formula 

5 = i?o(l + <Aw + <^s), 



(1) 



where Rq is the spread rate in the absence of wind, (j)^^ is the wind factor, and is the 
slope factor. The components of ([T]l computed from the fuel properties (Table [T}, the wind 
speed U, and the terrain slope tancp following the equations in Table |2] See [Rothermel 
( |1972[ ) for further details, derivation, and justifications. 
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Chaparral is a special fuel in that the spread rate depends only on wind speed. For 
chaparral, ([T]l is replaced by ( |Coen et aL||2001 Eq. 1) 



S = min{6,max{0.03333, 1.2974C/}} . (2) 



The only differences from Rothermel ( 1972 1 are the subtraction of the moisture from 



the fuel load in the computation rather than up front, limiting the slope and the windspeed, 
the special chaparral spread rate from CAWFE ([2]), and the explicit reduction of wind from 



6. 1 m height to midflame height, following Baughman and Albini ( 1980 1. 



In either case, the spread rate can be written as 

S = max|5o,i?o + cmin{e,max{0,C/}}'' + (imax{0,tan(/)}^| , (3) 

where ^o, Rq, b, c, d, e are the fuel-dependent coefficients that represent the spread rate 
internally. These coefficients are stored for every grid point. 

At a point on the fireline, denote by n the outside normal to the fire region, U the wind 
vector, and z the terrain height. The normal component of the wind vector, U = \J n, and 
the normal component of the terrain gradient, ta.n<p = Vz • n, are used to determine the 
spread rate, which is interpreted as the spread rate in the normal direction n. 

2.3 Fuel bumed and heat released 

Each location starts with fuel fraction F = 1. Once the fuel is ignited at a time ti, the fuel 
fraction decreases exponentially, 

F(t)=exp(-^^), t>U, (4) 

where t is the time, ti is the ignition time, Fq is the initial amount of fuel, and Tf is the 
fuel burn time, i.e., the number of seconds for the fuel to burn down to 1/e 0.3689 of the 
original quantity. Since by definition of the fuel weight w (Table[T]), the fuel burns down to 
0.6 of the original quantity in 600 s when w = 1000, we have 

(t-tj) 1000 / (t-ti) 

0.6 600 m =exp 

\ T{ 

which gives 

„ 600w w 



lOOOlnO.6 0.8514 



The input coefficient w is used in WRF-Fire rather than Tf for compatibility with existing 
fuel models and literature. 

The average sensible heat flux density released in time interval {t,t + At) is computed as 



, F{t)-F{t + At) 1 , -2 -n 

<^h = ^^ 7-rWeh, (Jm ^) (5) 

and the average latent heat (i.e., moisture) flux density is given by 

, F{t) -F{t + At) Mf + 0.56^ , _2 _n 

(^Q = ^^ -—r — ^ — Lwi, (Jm M (6) 

where 0.56 is the estimated mass ratio of the water output from the combustion to the dry 
fuel, and L = 2.5 x 10® Jkg~^ is the specific latent heat of condensation of water at °C, 
used for nominal conversion of moisture to heat. This computation is from CAWFE. 

3 Domain, grids, and nodes 

The atmospheric model operates on a logically quadrilateral 3-D grid on the Earth surface, 
and uses a sequence of horizontally nested grids, called domains ( Kalnay[ 2003| l. Only the 



innermost (the finest) atmospheric domain is coupled with the fire model; see also Sect. [8J 
Scalar variables in the atmospheric model are located at the centers of the 3D grid cells, 
while the wind vector components are at a staggered grid at the midpoints of the cell faces. 
The fire model operates on a refined fire mesh called the subgrid (Fig. [TJ, and all of its 
variables are all represented by their values at the centers of the cells of this fire subgrid. 

4 Mathematical core of the fire model 



Subsections 4.1 and 4.3 below follow Mandel et al. (20091. 



4.1 Fire propagation by the level-set method 

The model maintains a level-set function tp, the time of ignition t^, and the fuel fraction F. 
Denote a point on the surface by x = {x,y). The burning region at time t is represented by 
a level-set function ip = i]j{'x.,t) as the set of all points x such that ^(x,t) < 0. There is no 
fire at x if > 0. The fireline is the set of all points x such that = 0. Since on 
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the fireline, the tangential component of the gradient Vip is zero, the outside normal vector 
at the fireline is 



Now consider a point x(i) that moves with the fireline. Then the fire spread rate 5 at x 
in the direction of the normal n is 

S = n.|. (S) 

and, from the definition of the fireline, ?/;(x(t),t) = 0. By the chain rule and substituting 
from (|7]) and ([8]), we have 

So, the evolution of the level-set function is governed by the partial differential equation 
^ + S||VV||=0, (10) 



called the level-set equation (Osher and Fedkiw 2003 1. The spread rate S is evaluated from 
([3]) for all X, not just on the fireline. Since 5 > 0, the level-set function does not increase 
with time, and the fire area cannot decrease, which also helps with numerical stability by 
eliminating oscillations of the level-set function ifj in time. 

The level-set equation is discretized on a rectangular grid with spacing (Ax, Ay), called 
the fire grid. The level-set function ijj and the ignition time ti are represented by their values 
at the centers of the fire grid cells. This is consistent with the fuel data given in the center 
of each cell also. 

To advance the fire region in time, we use Heun's method (Runge-Kutta method of order 

2), 

= + At (^") + (v^"+i/2) ^ , (11) 

The right-hand side F is a discretization of the term — S || Vt/^H with up winding and artificial 
viscosity, 

F(V') = -5(U-n,Vz-n)||VV^|| +eAV', (12) 
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where n = VV'/H V'f/'ll is computed by finite central differences and Vijj = [V^V'; Vj^V'] 
is the upwinded finite difference approximation of Vip by Godunov's method ( |Qsher and 
Fedkiw||2003t p. 58), 







if V^?/'<0 and Vj'i/'<0, 
if V;;:V'>0 and V^^>0, 
if V^V'<0 and V^^>0, 



(13) 



otherwise V^. ip 


if 




> 






if 




< 





where V+V ^rid V^, are the right and left one-sided finite differences 

il){x + /\x,y)-il}{x,y) 



Ax 

'il;{x,y)-ip{x-Ax,y) 
Ax 



and similarly for Vyip and V^^ V- Further, in (12 1, e is scale-free artificial viscosity (e = 0.4 
here), and 



■ilj{x + Ax,y) — 2i{j{x,y) + ip {x — Ax ,y) 
~ Ax 



+ similar term for y 



is the five-point Laplacian of tp scaled so that the artificial viscosity is proportional to the 
mesh step, 

A^^Ax^ + Ay^. 



dx"^ 



dy^ 



A numerically stable scheme with upwinding, such as ( 13 1, is required to compute the 
term ||VV'|| in the level set equation (lOl. However, in our tests, the gradient by standard 
central differences. 



ip{x + Ax,y) — ip{x — Ax,y) ip{x,y + Ay) — ip{x,y — Ay) 



2Ax 



2Ay 



worked better in the computation of the normal vector n by (|7]), which is used to evaluate 
the normal component of the wind and the slope in Q. 



Before computing the finite differences up to the boundary, the level-set function is ex- 
trapolated to one layer of nodes beyond the boundary. However, the extrapolation is not 
allowed to decrease the value of the level-set function to less than the value at either of the 
points it is extrapolated from. For example, when is the last node in the domain in the 
direction x, the extrapolation 

tpi+ij = max{tpij + {^pij - tpi-ij) ,tpij,ipi-ij} , 

is used, and similarly in the other cases. This is needed to avoid numerical instabilities at 
the boundary. Otherwise, a decrease in V' at a boundary node, which may happen with non- 
homogeneous fuels in real data, is amplified by the extrapolation, and ip keeps decreasing 
at that boundary node in every time step until it becomes negative, starting a spurious fire. 

The model does not support fire crossing the boundary of the domain. When if; <0 is 
detected near the boundary, the simulation terminates. This is not a limitation in practice, 
because the fire should be well inside the domain anyway for a proper response of the 
atmosphere. 

4.2 Computation of the ignition time 

The ignition time ti in the strip that the fire has moved over in one time step is computed by 
linear interpolation from the level-set function. Suppose that the point x is not burning at 
time t but is burning at time t + At, that is, (x, t ) > and ip(^K,t + At) <0. The ignition 
time at x satisfies t/j (x,ti (x)) = 0. Approximating iphy a. linear function in time, we have 



V;(x,ti)-V;(x,t) _^(x,t + At)-V.(x,ti) 
ti{x)-t ^ t + At-tiix) 



and we take 



ti(x) = t4 



V'(x,t)At 



(14) 



ip{x,t)-i){:xi,t + At)' 



4.3 Computation of the fuel fraction 



The fuel fraction is approximated over each fire mesh cell C by integrating Q over the fire 
region. Hence, the fuel fraction remaining in cell C at time t is given by 




(15) 



xec 

i/'(x,t)<0 
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Once the fuel fraction is known, the heat fluxes are computed from Q and This scheme 
has the advantage that the total heat released in the atmosphere over time is exact, regardless 
of approximations in the computation of the integral ([T5]). Our objective in the numerical 



evaluation of ( [15] ) is a method that is second order accurate when the whole cell is on fire, 
exact when no part of the cell C is on fire (namely, returning the value one), and provides a 
natural transition between these two cases. Just like standard schemes in numerical analysis 
can be derived from the requirement that they are exact for all polynomials up to a given 
degree, the guiding principle here is that the scheme should be exact in as many special 
cases as possible. Then we expect that the scheme should work well overall. 

While the fuel burn time Tf can be interpolated as constant over the whole cell, the level- 
set function ip and the ignition time ti must be interpolated more accurately to allow a 
submesh representation of the burning area and a gradual release of the heat as the fireline 
moves over the cell. In addition, we need the fuel fraction computed over each mesh cell, 
because the heat fluxes in the mesh cells are summed up to give the heat flux in an atmo- 
spheric cell. Our solution is to split each cell into 4subcells Cj, interpolate to the corners 
of the subcells, and add the integrals, 

xec i— IxeCj 

i>{^)<0 i/'(x)<o 

cf., Fig.|2] The level-set function ip is interpolated bilinearly to the vertices of the subcells 
Cj, and the burn time Tf is constant on each Cj, given by its value at the fire grid nodes. 
However, to interpolate the ignition time ti we first define ti outside of the fire region and 
on the fireline by 

ti = t if V'>0. (17) 
This allows us to omit the condition < in the definition of the integration domains in 



(16 1 and integrate on the whole cells, respective subcells, only. Then, we interpolate ti 



bilinearly to the vertices of the subcells Cj and correct the resulting values by applying the 



compatibility condition ( 17 1. 



To compute the integral over a subcell Cj, we first estimate the fraction of the subceU 
that is burning, by 

areajx £ Cj : V'(x) < 0} ^ 1 / ELiV^(xfc) \ 

where are the the corners of the subcell Cj. This approximation is exact when no part 
of the subcell Cj, is on fire, that is, all if) (x^,) > and at least one (xjt) > 0; the whole Cj 
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is on fire, that is, all ipi'x.i^) < and at least one ip{^k) < 0; or the values 'ip{'Xk) define a 
linear function and the fireline crosses the subcell diagonally or it is aligned with one of the 
coordinate directions. 

Next, replace ti(xfc) by t when V'(xfc) > (i.e., the node is not on fire), and compute 
the approximate fraction of the fuel burned as 



(19) 



x6C 
i/.(x,t)<0 



This calculation is accurate asymptotically when the fuel bums slowly and the approxima- 
tion P of the burning area is exact. 

4.4 Ignition 

Typically, a fire starts from a horizontal extent much smaller than the fire mesh cell size, 
and both point and line ignition need to be supported. The previous ignition mechanism 



( Mandel et al.[ 2009 1 ignited everything within a given distance from the ignition line at 



once. This distance was required to be at least one or two mesh steps, so that the initial fire 



is visible on the fire mesh, and the fire propagation algorithm from Sect. 4.1 can catch on. 
This caused an unrealistically large initial heat flux and the fire started too fast. 

The current ignition scheme achieves submesh resolution and zero-size ignition. A small 
initial fire is superimposed on the regular propagation mechanism, which then takes over. 
Drip-torch ignition is implemented as a collection of short ignition segments that grows at 
one end every time step. Multiple ignition segments are also supported. 

The model is initialized with no fire by choosing the level-set function Tp (x,io) =const > 
0. Consider an initial fire that starts at time tg on a segment a,b and propagates in all 
directions with an initial spread rate until the distance rg is reached. At the beginning of 
every time step t such that 

we construct the level-set function of the initial fire, 

V'g(x,t) = dist(x,a7b) -5g(t-tg) (20) 
and replace the level-set function of the model by 



V'Cx,*) :=min{V'(x,f),^/;g(x,t)}. 
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(21) 



For a drip-torch ignition starting from point a at time tg at velocity v until time ih> the 
ignition line at time t is the segment a,a + v(min{t,th} — *g)> and (20l becomes 



V'g (x,t) = dist ( x,a,a + v(min{t,th} — *g) I ~™in{^g)'S'g(i — ig)} 



followed again by (21 1, at the beginning of every time step begining at time t such that 



trr<t<th + 



The ignition time of newly ignited nodes is set to the arrival time of the fire at the spread 
rate Sg from the nearest point on the ignition segment. 



5 Atmospheric model 



We summarize some background information about WRF-ARW from Skamarock et al. 



(2008 1, to the extent needed to understand the coupling with the fire module. 

The model is formulated in terms of the hydrostatic pressure vertical coordinate r], scaled 
and shifted so that r/ = 1 at the Earth surface and r/ = at the top of the domain. The 
governing equations are a system of partial differential equations of the form 

-^=Rm. (22) 

where R contains also the advection terms, and ^ = {U,V,W,(j)' ,@,iJ,' ,Qm)- The fundamen- 
tal WRF variables are ;U = fi{x,y), the hydrostatic component of the pressure differential 
of dry air between the surface and the top of the domain, written in perturbation form /i = 
]! + fi', where JI is a. reference value in hydrostatic balance; U = fiu, where u = u{x,y,r]) 
is the Cartesian component of the wind velocity in the x-direction, and similarly V and 
W; Q = 1x0, where 6 = 6{x,y,7]) is the potential temperature; (j) = (j){x ,y ,7]) = (j) + (j)' is the 
geopotential; and Qm = fJ-Qm is the moisture content of the air. The variables in the state 



<I> evolved by (22i are called prognostic variables. Other variables computed from them, 
such as the hydrostatic pressure p, the thermodynamic temperature T, and the height z, 
are called diagnostic variables. The variables that contain /i are called coupled. The value 



of the right-hand side R{^) is called tendency. See Skamarock et al. (2008 p. 7-13) for 
details and the form of R. 
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The system ( 22 1 is discretized in time by the explicit 3rd order Runge-Kutta method 
At 



$1 



3 

$2 = $* + ^i2($l) 



(23) 



where the differential operator R is discretized by finite differences and the tendencies 
from physics packages, such as the fire module, are updated only the third Runge-Kutta 



step (Skamarock et al. 2008 p. 16). In order to avoid small time steps, the tendency in the 
third Runge-Kutta step also includes the effect of substeps to integrate acoustic modes. 



6 Coupling of the fire and the atmospheric models 

The terrain gradient is computed from the terrain height at the best available resolution and 
interpolated to the fire mesh in preprocessing. Interpolating the height and then computing 
the gradient would cause jumps in the gradient, which affect fire propagation, unless high- 
order interpolation is used. 

In each time step of the atmospheric model, the fire module is called from the third step 
of the Runge-Kutta method. First the wind is interpolated to a given height zj above the 
terrain (currently, 6. 1 m following BEHAVE), assuming the logarithmic wind profile 

fconstlnl-, z>zo, 
\0 0<z<zo, 

where z is the height above the terrain and zq is the roughness height. For a fixed horizon- 
tal location, denote by zi, Z2,... the heights of the centers of the atmospheric mesh cells; 
these are computed from the geopotential (p, which is a part of the solution. The horizon- 
tal wind component u{zf ) under the n-points (Fig. [T]l is then found by vertical log-hnear 
interpolation, that is, u{zf) is found by 1-D piecewise linear interpolation of the values 
u{zq) = 0, u{zi), u{z2),... at Inzo, Inzi, lnz2,--- to Inzf. If Zf < zq, we set u{zf) = 0. 
The V component of the wind is interpolated vertically in the same way. Each horizontal 
wind component u, v is then interpolated separately to the cell centers of the fire subgrid 
by bilinear interpolation. 

The fire model then makes one time step: 

1. If there are any active ignitions, the level-set function is updated and the ignition times 
of any newly ignited nodes are set following Sect. |4.4| 
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2. The numerical scheme ( 1 1 1-( 13 1 for the level set equation ( 10 1 is advanced to the next 
time step. 

3. The time of ignition set for any any nodes that were ignited during the time step, from 



(14i 



4. The fuel fraction is updated following Sect. |4.3| 

5. The sensible and latent heat flux densities are computed from Q and Q in each fire 
model cell. 

6. The resulting heat flux densities are averaged over the fire cells that make up one at- 
mosphere model cell, and inserted into the atmospheric model, which then completes 
its own time step. 

The heat fluxes from the fire are inserted into the atmospheric model as forcing terms 
in the differential equations of the atmospheric model into a layer above the surface, with 
assumed exponential decay with altitude. Such scheme is needed because WRF does not 
support flux boundary conditions. This is code originally due to Clark et al. ( 1996a|b I and 



it was rewritten for WRF variables in Patton and Coen| ( |2004[ ). The sensible heat flux is 



inserted as an additional source term to the equation for the potential temperature 6, equal 
to the vertical divergence of the heat flux, 

— — (x,y,z) = i?0($) + ^ — exp , 

at ag{x,y,z) dz \ Zcxt/ 



where Rq{^) is the component of the tendency in the atmospheric model equations (22), 
a is the specific heat of the air, g{x,y,z) is the density, and Zext is the heat extinction 
depth, given as parameter f ire_ext_grnd in namelist . input. The latent heat flux 
is inserted similarly into the tendency of the vapor concentration q^^ by 

— i.,y,z) = R^^ m + ^^(,^^^,) ^exp [-—) , 



where L is the specific latent heat of the air. Cf. Clark et al. ( 1996a[ Eqs. 10, 12, 13, 18). 



7 Software structure 
7.1 Parallel structure 

Parallel computing imposes a significant constraint on user programming technique. WRF 
parallel infracture ( ,MichalakeS| |2000) divides the domain horizontally into patches. Each 
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patch executes in a separate MPI process and it may be further divided into tiles, which 
execute in separate OpenMP threads (Fig. |3]l. Communication between the tiles is accom- 
plished by exiting the OpenMP parallel loop over the tiles. The fire grid tiles are colocated 
with the atmospheric grid tiles (Fig. [T]). The patches are declared in memory with larger 
bounds than the patch size, and communication between the patches is accomplished by 
HALO calls (actually, includes of generated code), which update a layer of array entries 
beyond the patch boundary from other patches. The fire module computational code itself 



is designed to be tile-callable as required by the WRF coding conventions (WRF Working 



Group 2[ 2007 1. Tile-callable code updates array values on a single tile, assuming that it 



can safely read data from a layer of several array entries beyond the tile boundary. The 
communication (OpenMP loops or HALO calls) happens outside; this means that every 
time when communication is needed, tile-callable code must exit, and then the next stage 
can resume on the next call (Fig. |4]l. The fire module code executes in 6 stages interleaved 
with communication, 3 stages for initialization and 3 stages in every time step. 

7.2 Software layers 

The fire module software is organized in several isolated layers (Fig. [5]). The driver layer 
contains all exchange of data between the tiles in parallel execution. The rest of the code 
is tile-callable. The driver layer calls the interpolation and other coupling between the fire 
and the atmospheric grids, and the fire code itself. The atmospheric physics layer mediates 
the insertion of the fire fluxes into the atmosphere, as described in Sect. [2] Only these 
two layers depend on WRF; the rest of the fire module can be used as a standalone code, 
independent on WRF. The utility layer contains interpolation and other service code, such 
as stubs to control access to WRF infrastructure, so that WRF calls can be easily emulated 
in the standalone code. The model layer is the entry point to the fire module. The core 
layer is the engine of the fire model, described in Sect. [4] The fire physics layer evaluates 
the fire spread rate and heat fluxes from fuel properties. One of the goals of the design is 
that the only components that will need to be modified when the fire module is connected 
to another atmospheric model in future are the driver layer, the atmospheric physics layer, 
and the WRF stubs in the utility layer. 
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8 Recommended WRF settings 
8.1 Domains and nesting 

WRF-Fire may be run in both "ideal" and "real" modes, which require slightly differ- 
ent setups. In both cases, the model requires a set of data defining model initialization 
(wrf input). In the real cases, boundary conditions in a form of wrfbdy files must be 
also provided, and both types of files are created by real . exe preprocessor from the WRF 
Preprocessing System (WPS). These files contain not only meteorological and topograph- 
ical data but also fire related information, such as the fuel type map and high-resolution 
topography on the fire mesh. Since the WRF-Fire initialization for the real cases does not 
differ from the one for the regular WRF, all physical and dynamical options available in the 
regular WRF are also available in WRF-Fire. Therefore, the same general rules apply to the 
configuration of WRF-Fire as to the configuration of the regular WRF. However, one should 
keep in mind that resolutions of the finest domains in fire simulations are usually signifi- 
cantly higher than in weather forecasting applications. This has two consequences in terms 
of the proper WRF-Fire setup. First, if the resolution of any of the inner domains is less than 
100 m, this domain should be actually resolved in the large eddy simulation (LES) mode, 
without the boundary layer parameterizations. At this resolution, the model should be able 
to resolve the most energetic eddies responsible for mixing within the boundary layer, so 
the boundary layer parameterization in this case is not needed. Second, since in the nested 
mode, vertical levels are common for all domains, the height of the first model level se- 
lected for the most outer (parent) domain, defines also the level of the first model layer 
for all inner (child) domains, even if their horizontal resolutions are an order of magnitude 
smaller. The fact that the vertical model resolution is the same for all domains significantly 
limits the minimum height above the ground of the first model level. This in turn is cru- 
cial for the fire model, which uses the wind speed interpolated to 6. 1 m above the ground. 
Therefore, in the cases when the first model level must be relatively high above the ground 
it is recommended to perform downscaling using the ndown .exe program, being a part 
of the WRF distribution. In this case the outer domains are run separately without the fire, 
and then based on the output from this simulation, ndown . exe creates a set of new initial 
and boundary condition files (wrf input and wrfbdy) for the separate simulation from 
the innermost domain(s). This allows for a new setup of vertical levels for the innermost 
domains, and selecting proper physical options for them. 
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8.2 Large Eddy Simulation and surface properties 

To enable the high-resolution simulation in Large Eddy Simulation (LES) mode, user should 
first disable the boundary layer parameterization (bl_pbl_physics=0). The LES mode 
requires the proper surface fluxes in order work properly. We recommend the option 
isf f lx=l, which makes WRF use a surface model to compute the surface fluxes. Other 
options with constant heat fluxes and drag are not well suited for fire simulations. Out of all 
surface exchange parameterizations only the classic Monin-Obukhov theory (s f _s f c 1 ay_phy s i 
is recommended for the LES cases. This option assures a proper computation of surface 
transfer coefficients that are used together with the surface properties (provided by the sur- 
face model) for computation of the surface fluxes of the momentum, heat and moisture. The 
surface model itself computes properties of the surface, but does not compute the surface 
exchange coefficients, which are needed for computation of the surface fluxes. Hence, in 
order to compute them, the surface properties must be provided by a surface model, which 
is enabled by choosing a non-zero sf _surf ace_physics. The subgrid scale parameter- 
ization used by the WRF in LES mode is defined by the km_opt parameter, which should 
be set to 2 (TKE closure), or 1 (Smagorinsky scheme). 

In real cases, real . exe automatically provides proper initialization for the selected 
land surface model, and all other components. In idealized cases, users have an option 
of the basic surface initialization, intended to be used without the surface model, or the 
full surface initialization (sf c_f ull_init=l). One should keep in mind that without 
the full surface initialization, there is no direct way to define surface properties such as 
temperature or roughness. For idealized cases with the full surface initialization, the surface 
scheme utilizes a table containing records of land-use categories and corresponding surface 
properties like roughness length, heat capacity, etc. All these properties are defined in a text 
file LANDUSE . TBL, which may be edited by the user. Therefore, setting up the land-use 
category is enough to provide all static surface properties. The basic parameters required by 
this surface model like land use index, surface air temperature and soil temperature, may be 
defined directly in namelist . input by the variables sf c_lu_index, sf c_tsk, and 
s f c_tmn if they are intended to be the same over the whole domain. If they are not spatially 
uniform, they may be read in from external files if f ire_read_lu, f ire_read_tsk, or 
f ire_read_tmn are set to true. For details about the input data for real cases, see Sect.|9j 

8.3 Fire subgrid refinement ratios 

The fire mesh needs to be about 10 times finer than the atmospheric mesh to allow for 
gradual heat release into the atmosphere, even if fuel and topography data may not be 
available at such fine resolution. The fire mesh refinement in the x and y direction (sr_x 
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and sr_y) must be defined in the domain section of namelist . input. Since these 
refinement factors define dimensions of the fire-related variables, they must be selected 
before execution of real.exe, which generates the WRF input files. Any change in 
atmospheric to fire grid ratios requires re-running real . exe and creating new input files. 
The atmospheric mesh step should be about 60 m or less for proper feedback of the wind 
on the fire line; larger mesh step can result in incorrect fire spread rates and atmospheric 
behavior ( |Clark et al.|[T996a| p. 887). 



8.4 Time step 

In real WRF-Fire simulations performed in multi-domain configurations the time step re- 
quirements for the outer domains (run without fire) do not differ from general meteoro- 
logical cases. The recommended time step of 6 times the horizontal grid spacing (in km) 
may be used as a starting point. However, for the finest domains run with fire simulations, 
the time step in most cases must be significantly smaller. For domains with low vertical 
resolution and simple topography, the horizontal mesh step is crucial for numerical stabil- 
ity, since the horizontal velocity is greater than the vertical one. In fire simulations with 
high vertical resolution, the vertical velocity induced by fire may violate the CFL condi- 
tion. Therefore, it is advisable to use a vertically stretched grid, with finer resolution at 
the surface (where updraft velocities are not that high) and lower resolution at higher levels 
where stronger updrafts are expected. This allows for having the first model level relatively 
close to the ground, yet with vertical spacing aloft big enough to handle strong convective 
updrafts without violating the CFL condition. 

In real cases, the pressure levels may be defined directly in the namelist . input 
file. In ideal WRF-Fire runs, there is now an option stretch_hyp, which turns on hy- 
perbolic grid stretching. The grid refinement may be adjusted using the z_grd_scale 
namelist variable. One should keep in mind that running the WRF-Fire simulations with 
high-resolution topography in most cases limits the maximum numerically stable time step. 
Steep terrain often induces high vertical velocities that may violate the CFL condition. 
Therefore, these cases usually require significantly smaller time steps than similar simula- 
tions run with low-resolution, smooth topography. 

9 Data input 

A WRF ideal run is used for simulations on artificial data. An additional executable, 
ideal . exe, is run first to create the WRF input. A different ideal . exe is built for 
each ideal case, and the user is expected to modify the source of such ideal case to run 
custom experiments. The ideal run for fire supports optional input of gridded arrays for 
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land properties, such as terrain height, roughness height, and terrain height. This allows 
to run simulations which go beyond what would normally be considered an ideal run and 



simplifies custom data input; the simulation of the FireFlux experiment (Sect. 10) was done 
in this way. 

A WRF real run is used for prediction and analysis of natural events. For a real run, a 
user must supply data for the initial and boundary conditions for the WRF simulation. The 
WRF Preprocessing System (WPS) ( |Wang et al.[ [2010 , Chapter 3) contains a number of 
utilities useful for preparing standard atmospheric and surface datasets for input into WRF. 

1. Geogrid creates the surface mesh from a specified geographic projection and interpo- 
lates static surface data onto the mesh. It supports several interpolation methods as 
well as data smoothing and creation of gradient fields. Geogrid reads data in a tiled 
binary format described by a text file and writes to a NetCDF file for each nested 
mesh. All data required for atmospheric simulations up to 30 arcseconds resolution 
globally are provided by NCAR. 

2. Ungrib extracts atmospheric data from standard GRIB files and writes to a simple bi- 
nary format. Ungrib does not do any interpolation; it only searches through a number 
of files for necessary variables within the time window of the simulation. Data for 
ungrib must be obtained by the user. Several free sources of atmospheric GRIB data 
are available online from production weather simulation. 

3. Metgrid reads the output from geogrid and ungrib and produces a series of NetCDF 
files read by WRF's real . exe binary. The geogrid output is copied directly into 
each of these files, while the ungrib output is interpolated horizontally on to the com- 
putational mesh. 

The metgrid files produced by WPS are portable and relatively compact so they can 
be transferred to a computer cluster for the simulation's execution. From this point, the 
real . exe program in WRF handles the vertical interpolation of atmospheric fields and 
all processing for the creation of WRF's initial (wrf input) and boundary (wr f bdy) files. 

WPS has been extended with the ability to produce data defined on the refined surface 
meshes used by WRF-Fire (Sect. [8]l; however, it is not possible to distribute high resolu- 
tion, global fields as is done in the standard dataset. Instead, the user must download any 
necessary high resolution fields and convert them into geogrid's binary format for each sim- 
ulation. WRF-Fire is distributed with an additional utility, convert_geotif f . x, which 
can perform this conversion from any GeoTIFF file. This utility is written entirely in C and 
depends only on the GeoTIFF library. 

For a WRF-Fire simulation, it is only strictly necessary to download one additional 
dataset for input into geogrid. This dataset contains fuel behavior categories and is stored 
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in the variable NFUEL.CAT. For simulations within the United States, this data can be ob- 
tained in GeoTIFF format from the USGS at http://www.landfire.gov. WRF-Fire uses an 
additional variable for topography, ZSF, which is allowed to be different from the topog- 
raphy used used by the atmospheric code defined by HGT. This is useful because a high 
resolution WRF simulation generally requires the topography to be highly smoothed in 
preprocessing for numerical stability. The fire code can benefit from a rougher topography 
for more accurate fire spread computations. 

Once the static data is converted into the geogrid binary format, the GEOGRID . TBL 
should be edited to inform geogrid of the location of each supplementary dataset. WRF-Fire 
expects two variables to be created on the refined subgrid (NFUEL_CAT and ZSF), this is 
indicated by the line subgrid=Yes; all other variables will be defined on the standard 
atmospheric grid. 

For atmospheric data, it is best to use the highest resolution dataset available to initialize 
a WRF-Fire simulation to capture as much of the local conditions near the fire as possible. 
Generally, publicly available atmospheric data is limited to around 10 km resolution. As a 
consequence, one should create several nested grids, each with a 3 to 1 refinement ratio, 
and a long spin-up prior to ignition in order to recreate local conditions. Preliminary re- 
sults indicate that assimilation of data from weather stations or satellite radiances may be 



required for an accurate simulation ( ,Beezley et al. 20101. 



10 Computational simulations 



Kim ( 2011[ l has verified that the level-set method in the fire module advects the fire shape 
correctly, on some of the same examples that were used to verify the tracer code in CAWFE 
( |Clarketal.||2004l l. 

A number of successful simulations with WRF-Fire now exist. 

Jenkins et al.l(|2010]) have demonstrated fireline fingering behavior for a sufficiently long 



fireline (Figs. [8} |9]l on an ideal example, with similar results as in Clark et al. ( 1996a|b ) 



Kochanski et al.| ( j2010 1 have demonstrated the validity of WRF-Fire on a simulation of 
the Clements et al.| ( |2007 1 FireFlux grass fire experiment and obtained good agreement 
with data (Figs. |6j [7]). pobrinkova et al. ( |2010 l simulated a fire in Bulgarian mountains 



using real meteorological and geographical data, and ideal fuel data. Beezley et al. (20101 
simulated the 2010 Meadow Creek fire in Colorado mountains using real data from online 
sources. Topography (Fig. [TO]) at up to 3 m horizontal resolution was obtained from the 



National Elevation Dataset (NED, http://ned.usgs.gov ) and fire fuel datasets from Landfire 
(jhttp://landfire.cr.usgs.govll at up to 10 m resolution. Six nested domains were required 
to scale the simulation down from the atmospheric initialization (32 km) to the fire grid 
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resolution (10 m). Cloud physics was enabled in domains 1-3. The fire subgrid refinement 
ratio was 10 times on the finest domain to capture fire surface variables and for a gradual 
release of the heat flux near the fireline. Realistic fire and atmosphere behavior was obtained 
(Figs. 
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11 Discussion 

11.1 Additional features 

WRF-Fire does not yet support canopy fire, although canopy fire colocated with ground fire 
is contained in CAWFE. The reason was the desire to keep the code as simple as possible 
early on and add features only as they can be verified and validated. The support for canopy 
fire will be added in future. Adding smoke from the fire to WRF is also under consideration. 
A list of desired features and a record of the progress of the development are maintained at 



http://www.openwfm.org/wiki/WRF-Fire_development_notesj 
11.2 Atmospliere 

Rothermel's spread model ([T]l assumes wind as if the fire was not there. In practice, the 
wind was measured away from the fire. In a coupled model, however, the feedback on 
the fire is from the wind that is influenced by the fire. [Clark et aL ( 2004 1 noted that the 
horizontal wind right above the fireline may even be zero, and proposed to take the wind 
from a specified distance behind the fireline. Also, the strong heat flux from fire disturbs 
the logarithmic wind profile, and the rate of spread as a function of wind at a specific 
altitude may not be a good approximation; rather, the fire spread may depend more strongly 
on the complete wind profile ([Jenkins et al.[ [2010) and on turbulence (Sun et al.[ [2009[). 



The assumption of horizontal homogeneity in the Monin-Obukhov similarity theory is not 
satisfied here; the horizontal dimension of the active part of fire is not orders of magnitude 
larger than the boundary layer height as required, and it may be in fact smaller. Another 
indication that the Monin-Obukhov theory may not apply for fires is a strong drop in the 
heat transfer in the case of strong temperature gradients, shown in our preliminary tests. 

Horizontal wind could be interpolated vertically to different heights for different fuels 
like in CAWFE model, which takes the wind from different mesh levels for different fuels. 



However, here we follow a classical approach of Rothermel (1972i and Baughman and 



Albini (1980 1, where the wind speed is evaluated at the common 6.1m height, and then 



converted to the mid-flame height using the fuel-specific wind correction factors. 

Very strong vertical components of the wind caused by the fire result in the need for short 



time steps to avoid violation of the vertical CFL condition (Sect. 8.4 1. It would be interest- 
ing to couple the fire module also with the Non-hydrostatic Mesoscale Model (NMM) core 
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of WRF, which is implicit in the vertical direction ( Janjic et"aL| 2005 1, and it may perform 
better in the presence of strong convection ( Litta and Mohanty] 2008 1. The ARW core is 
semi-implicit in the vertical direction in the vertical wind component and the geopotential. 

11.3 Fire 



The more recent Scott and Burgan| ( |2005| l fuel categories are more detailed than Anderson 
( |1982[ ) categories, they are supported by BehavePlus, and fuel maps using them are avail- 
able from Landfire. But instead of describing additional categories in namelist . fire, it 
may be more useful to support the import of fuel files from BehavePlus, which is also well 



suited for editing and diagnosing fuel models. More accurate fuel models (Albini et al. 



1995 ; Clark et al. 1996a l, including those in BehavePlus, consider fuels to be mixtures of 



components with different burn times, which results in a different heat release curve. 

While the spread rate of established fire in the simulation of the FireFlux experiment 
was reasonably close, the simulated fire still arrived at the observation towers too soon 



( [Kochanski et al.[ |2010| ), because it started too quickly. A better parametrization of the 
ignition process seems to be in order. The fire spread in the Meadow Creek fire simulation 
was also too fast, but for a different reason. It is well known that the actual spread rates 
of wildland fires tend to be lower than the spread rates in simulations, which are derived 
from laboratory experiments. This effect might be attributed to irregularities on scales not 
captured by the simulation (Finneyj 1998 p. 34), including granularity of the fuel supply 
not reflected in the data. Refining the semi-empirical model from detailed numerical simu- 
lations and parametrizing complex fire behavior are suggested important research areas. 

The computation of the heat fluxes in (|5]l and Q does not take into account the evap- 
oration of moisture present in the fuel, only the production of water by burning of hydro- 
carbons. This error is typically just few %, however, which is small in comparison with 
other uncertainties. The fuel models should be dynamic (with variable fuel moisture) as in 



BehavePlus. Coen (20051 added an explicit diurnal cycle for the moisture into CAWFE. 
Here, moisture content could be coupled with existing WRF land surface models, which 
could take into account air humidity and precipitation. The radiative and convective parts 
of the sensible heat flux should be treated differently. The release of surface heat and mois- 
ture into the atmosphere are already present in WRF soil models. Their scale, however, is 
different from the powerful heat release from a fire. 

11.4 Numerical methods 

In a numerical implementation, the level-set method is global, unlike tracers, which move 
locally. In spite of the fact that the level-set equation determines the fire spread locally from 
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the spread rate at the fireUne, the behavior of the fireUne depends shghtly on the wind, the 
fuel, and the level set function in certain other locations from previous time steps, because 
of the discretization errors and the artificial diffusion. This nonlocal behavior has not been 
practically significant, however. 



The fuel fraction calculation ( 19 1 can have significant error in the fire subgrid cells near 
the fireline, which will to some degree average out over the atmospheric mesh cells. Rig- 
orous error analysis will be done elsewhere. We are currently testing an alternative method 
which is always second order in the sense that it is exact when the time from ignition and 
the level-set function are linear in space. The alternative method is more computationally 
expensive, but, on the other hand, it might allow to decrease the subgrid refinement ratio; 
with large meshes, it is possible to run against 32 bit integer limits. 

11.5 Data assimilation 



Data assimilation for wildland fires is an area of great interest. Methodologies for a reaction- 
diffusion model were proposed based on the ensemble Kalman filter (EnKF) and the particle 



filter ( |Mandel et al.[|2004i ). Unfortunately, statistical perturbations can cause spurious fires, 
which do not dissipate. Combination of the EnKF with Tikhonov regularization alleviates 



the problem somewhat (Johns and Mandel 2008 Mandel et al.[ 20091, but the resulting 
method is still not robust enough. A new method, called morphing EnKF and based on 
combined amplitude and displacement correction ( Beezley and Mandel} 2008"), was shown 
to work with WRF-Fire ( Mandel et al.[ 2009 1, and it is under continued development (Man^ 
del et al.[ 2010[ 2011[ ). We are not aware of any work elsewhere on data assimilation for 
a coupled fire-atmosphere model. Particle filters were proposed for discrete cell-based fire 
models CBianchini e t al.[ 2006 [ |Gu et al.[ 2009 1, using fitness functions involving the area 
burned rather than intensities of physical variables. 

Starting the model from a known fire perimeter is important for many potential users. 
This can be understood as a data assimilation problem, but we are considering a simpler 
method for this particular case: prescribe the fire history up to the time of the given perime- 
ter to allow the atmospheric conditions to evolve, then allow the coupled model take over. 
Tools to produce such artificial fire history are being developed. Possibly the simplest al- 
ternative is an interpolation from a given ignition point and time to the given perimeter. A 
more complex version would run the fire model (without atmosphere) backwards in time 
and attempt to find the ignition point automatically. The latter approach could be also in- 
teresting for forensic purposes. 
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12 Conclusions 



We have described the coupled atmosphere-fire model WRF-Fire. The software is publicly 
available and it supports both ideal and real runs. Visualization and diagnostic utilities are 
available. Currently, the model is suitable for research and education purposes. Validation 
is in progress. 
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Table 1. Fuel properties. The notation is from'Rothermel ( 1972 1 except as indicated. The identifiers 
are as used in WRF-Fire and CAWFE. In the input files, some quantities are given in English units 



per |RothermeIl ( [T972l i; see | Wang etaL] ( |20T0{ p. A-5) 



symbol description 



identifier 



w 

a 

P-P 
St 

h 



wind adjustment factor ( Baughman and Albini 1980 1 

from 6. 1 m to midflame length 

fuel weight (i.e., burn time) (s) 

40% decrease of fuel in lOmin for w = 1000 

total fuel load (kg m~^) 

fuel depth (m) 

fuel particle surface-area-to-volume ratio (1/m) 

moisture content of extinction (1) 

ovendry fuel particle density (kg m^'^) 

fuel particle total mineral content (1) 

fuel particle effective mineral content ( 1 ) 

fuel heat contents of dry fuel (J kg^^) 

fuel particle moisture content (1) 



windrf 



weight 
f gi 

f ueldepthm 
savr 
f uelmce 
f ueldens 

St 

se 

cmbcnst 

f uelmc_g 
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Table 2. Computation of the fire spread rate factors in (T| from the fuel properties (Table[TJ, the wind speed U at 6.1 m, 
and the terrain slope tant/i. All equations are from Rothermel 1 19721 unless otherwise indicated. All input quantities are 
first converted from metric to English units (BTU-lb-ft-min) to avoid changing the numerous constants in the Rothermel] 
{1972} computations. Further, following CAWFE, the wind is limited to between and 30 ms ^ and the slope is limited to 
nonnegative values. 



equation 

exp[(0. 792+0. 681ct'''')(^+0.1)] 
^ ^ 192+0. 2595(T 

Ill = TwnhrjMVs 
r,,=0.174S'-0-i9 

r,M = l-2.59||i+5.11 



-3.52(1^ 



1+St 



max 495+0.5940-1 



r 

^ ma 

Pb = 

4.77O-0 I-7.27 

e = exp(-lf) 

Q,g = 250/3 + lllGMy: 

<^W = Cmax[/f (^)^ 
C = 7.47exp(-0.133(70-S5) 

Ua = aU 

E = 0.715exp (-3.59 x IQ-^a) 
(7is = 5.275/3"*'-^tanV 



description 


source 


spread rate without wind 


Eq. 


(52) 


propagating flux ratio 


Eq. 


(42) 


reaction intensity 


Eq. 


(52) 


mineral damping coefficient 


Eq. 


(30) 


moisture damping coefficient 


Eq. 


(29) 


fuel loading net of minerals 


Eq. 


(24) 


total fuel load net of moisture 


from CAWFE 


optimum reaction velocity 


Eq. 


(36) 


maximum reaction velocity, 


Eq. 


(36) 


packing ratio 


Eq. 


(31) 


oven dry bulk density 


Eq. 


(40) 




Eq. 


(39) 


effective heating number 


Eq. 


(14) 


heat of preignition 


Eq. 


(12) 


wind factor 


Eq. 


(47) 




Eq. 


(48) 



adjustment to midflame height 
slope factor 



Table [Tlhere 
Eq. M 
Eq. (51 Iq) 
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Fig. 2. Division of fire mesh cells into subcells for fuel fraction computation. The level-set function 
tp and the ignition time ti are given at the centers ai,...,a4 of the cells of the fire grid. The inte- 
gral ( 16 1 over the cell C with the center is computed as the sum of integrals over the subcells 
Ci,...,C4. While the values of ijj and ti are known at =X3, they need to be interpolated to the 
remaining comers xi, X2, X4 of the subcell Ci from their values at the points ai , . . . ,34. 
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Fig. 3. Parallel communication in WRF. The computational domain is divided into disjoint rectan- 
gular patches. Each patch is updated by a single MPI process (distributed memory parallelism), and 
the process may read arary data in a strip around the patch, called halo region. The communication 
between the patches is by halo calls to the RSL parallel infrastructure ( Michalakes 2000[ l, which 



update the halo regions by the values from the neighboring patches. Each patch may be divided into 
tiles, which execute in separate OpenMP threads (shared memory parallelism). Following WRF 
coding conventions ( |WRF Working Group 2| |2007| l, computational kernels execute in a single tile. 
They may read array values from a strip beyond the tile boundary but no explicit communication is 
allowed. 3-D arrays are divided into patches and tiles in the horizontal plane, cf., Fig.[T] 
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WRF 



WRF 



Loop over fire 
code stages 



f 



* 4 

Halo exchange (MPI)<^JzZj^ Halo exchange (MPl) 
SOMP PARALLEL DO $OMP PARALLEL DO 

4 \ 4 \ 



Loop over fire 
code stages 



Fire model 
tile code 



Fire model 
tile code 



Fire model Fire model 
tile code tile code 



% 4 % 4 

$OMP END PARALLEL DO $OMP END PARALLEL DO 
Memory synchronized Memory synchronized 



WRF 



WRF 



Fig. 4. Parallel structure of the fire module in the WRF physics layer The core code itself executes 
on a single tile, with all communication done outside. Multiple passes through the fire module are 
needed in each time step. 
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WRF: Runge-Kuttastep3 



Driver: get grid variables, get 
flags, interpolation calls, 
OpenMP loops, OM halos 



Model: one time step, winds 
in, heat fluxes out 



Core: Advance the level set 
equation, ignition, compute 
fuel fraction, Dimensionless. 



Atmospheric physics: get 

temperature and moisture 
tendencies from heat fluxes 



Fire physics: sensible and 
latent heat fluxes from fuel 
fraction, fire rate of spread 



Utilities: interpolation, WRF stubs, debug I/O,... 
WRF: error messages, log messages, constants,... 



Fig. 5. Software layers of WRF-Fire. All physics dependencies are in the dashed box. The utilities 
layer is called from all the other layers above. 
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Fig. 6. Simulation of the FireFlux experiment (Clements et al. 20071 by WRF-Fire. Left: map of 



landuse category for the experimental plot, with the ignition line and the observation towers marked. 
Right: simulated and measured temperature profiles at the location of the observation towers. The 
simulated fire propagation takes 243 s from tower MT to tower ST, while the measured time is 255 s 
(4.7% difference). From ,Kochanski etaL] ( |20Tol i. 
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Fig. 7. Simulation of the FireFlux experiment ( Clements et al.| [2007) by WRF-Fire. Left: surface 
heat flux and selected flowlines. Visualization in VAPOR by Bedfich Sousedik. Surface image from 
Google Earth. Right: vertical velocity at 2 m height at tower S T. (See Fig.[6|left for location.) The 
simulation shows a good agreement with the experiment. From Kochanski et al. ( 2010) l. 
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Fig. 8. X-Y section of wind vector at 18 m and pressure perturbation 240 s after line ignition, 
initialized with uniform wind profile. The fire develops two fingers due to wind direction inversion 
in the middle. From [Jenkins et alT] ( |2010| l. 
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P_pert(kPa) Mir, UqK^(-0.36, 0.61x10"') 




□ 12 3 4 

Fig. 9. X-Z section of wind vector and pressure perturbation at the centerline for the fire in Fig.jSj 
From [Jenkins et al] ( |20T0| l . 
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Fig. 11. The finest domain in the Meadow Creek fire simulation 5h after ignition. Unburned fuel 
is displayed as green, burned fuel as brown. The heat flux from the fire appears near the fire line. 
Arrows indicate the surface winds, while streamlines show the atmospheric winds flowing over the 
fire region. Visualization in MayaVi. From Beezley et al. ( 2010[ l. 
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Fig. 12. The top level domain in the Meadow Creek fire simulation 5 h after ignition. Streamlines 
show the winds blowing East, over the Rocky Mountains and South down the coast of California. 
Visualization in MayaVi. From Beezley et al.| ( |2010| . 
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